use rand::distributions::Uniform;
extern crate nalgebra as na;
use na::*;
use nalgebra::{DMatrix, DVector};
use rand::prelude::Distribution;
use rand::seq::SliceRandom;
use rand::thread_rng;
use rand::Rng; // 0.8.5
use rand_distr::{Normal, NormalError};
use statrs::statistics::{MeanN, VarianceN};
use std::f64::consts::PI;
/*
先生成一个0到N间的随机整数i，代表选择第i列；

再生成一个0到1间的随机数，若其小于事件i占第i列矩形的面积的比例，则表示接受事件i，否则，接收第i列中不是事件i的另一个事件。

其实你可以算下这种方式事件i的概率，完全对应原来的概率分布。
*/

/// # Example
///
/// ```
/// use rspp::tools::resample;
/// let data = vec![1, 2, 3, 4, 5];
/// let weights = vec![0.4, 0.3, 0.2, 0.2, 0.6];
/// let ll = resample::resample_dup::<i32>(&data, &weights, 10);
/// dbg!("重采样===",ll);
/// ```
pub fn resample_dup<T: Clone>(data: &Vec<T>, weights: &Vec<f64>, size: usize) -> Vec<T> {
    let mut rng = rand::thread_rng();
    let n = data.len();
    let uni = Uniform::new(0.0, 1.0);
    let total_weight: f64 = weights.iter().sum();
    let norm_weight: Vec<f64> = weights.iter().map(|v| v / total_weight).collect();

    let mut res = Vec::new();
    while res.len() < size {
        let index = rng.gen_range(0..n);
        let threshold = uni.sample(&mut rng);
        if norm_weight[index] > threshold {
            let el = &data[index];
            res.push(el.clone());
        }
    }
    res
}
/// # Example
///
/// ```
/// use rspp::tools::resample::jacbin_one_row;
/// fn test_jacbin_f(x: &Vec<f64>) -> f64 {
///   x.iter().map(|i| i.powf(2.0)).sum()
/// }
/// let x: Vec<f64> = vec![1.0, 2.0];
/// let jac = jacbin_one_row(&x, test_jacbin_f);
/// dbg!("雅可比", jac);
/// ```
pub fn jacbin_one_row(x: &Vec<f64>, f: fn(&Vec<f64>) -> f64) -> Vec<f64> {
    let dt = 1e-8;
    let mut jac = Vec::new();
    for col in 0..x.len() {
        let djx = if x[col] != 0.0 { x[col] * dt.abs() } else { dt };
        let x1: Vec<f64> = x
            .iter()
            .enumerate()
            // 其他稳着不动,在col方向步进
            .map(|(i, v)| if i != col { *v } else { *v + djx })
            .collect();

        let res = (f(&x1) - f(x)) / djx;
        jac.push(res);
    }
    jac
}

pub fn random_from_stdnorm() -> f64 {
    let mut rng = rand::thread_rng();
    let normal = Normal::new(0.0, 1.0).unwrap();
    normal.sample(&mut rng)
}

pub fn gaussian_likehood(x: f64, std: f64) -> f64 {
    let n1 = -1. * x.powf(2.) / (2. * std.powf(2.));
    let n2 = (2. * PI).sqrt() * std;

    n1.exp() / n2
}

pub fn resample_cumsum(w: &Vec<f64>, np: usize) -> Vec<usize> {
    let wsum: f64 = w.iter().sum();
    let norm_w: Vec<f64> = w.iter().map(|v| v / wsum).collect();
    let w_cum: Vec<f64> = norm_w
        .into_iter()
        .scan(0.0, |acc, x| {
            *acc += x;
            Some(*acc)
        })
        .collect();
    let step = 1. / np as f64;
    let base: Vec<f64> = (0..np)
        .into_iter()
        .scan(0. - step, |acc, x| {
            *acc += step;
            Some(*acc)
        })
        .collect();

    let mut rng = rand::thread_rng();
    let uni = Uniform::new(0.0, step);
    let re_sample_id: Vec<f64> = base.iter().map(|v| v + uni.sample(&mut rng)).collect();
    let mut ind = 0;
    let mut index = Vec::new();
    for ip in 0..np {
        while re_sample_id[ip] > w_cum[ind] {
            ind += 1
        }
        index.push(ind);
    }
    index
}

pub fn rand_n_index(tlen: usize, n: usize) -> Vec<usize> {
    let l = 0..tlen;
    let mut numbers: Vec<usize> = l.into_iter().collect();
    numbers.shuffle(&mut rand::thread_rng());
    let res = numbers[0..n].to_vec();
    res
}
#[test]
fn test_resample() {
    let p = vec![1., 2., 3., 4.];
    resample_cumsum(&p, 100);
    let res = rand_n_index(6, 3);
    dbg!(res);
}
